# 0. initial setup ----
# clear workspace 
rm(list = ls()) 

# set seed
set.seed(11081984)

# set working directory
setwd("YOUR WORKING DIRECTORY")

# load required packages
library(tidyverse)
library(gridExtra)

# define background color
bg_colour <- "white"


# 1. load data frames with model coefficients ----
df <- read_csv("tables and figures/model_coef.csv") %>%
  filter(!var %in% c("cut1", "cut2", "cut3", "var(_cons[country_id])",
                     "0b.female", "0b.married")) %>%
  mutate(lower_95 = beta - 1.96 * se,
         upper_95 = beta + 1.96 * se, 
         lower_99 = beta - 2.58 * se,
         upper_99 = beta + 2.58 * se,
         var = recode(var,
                      "1.religion" = "none",
                      "2b.religion" = "roman catholic (ref.)",
                      "3.religion" = "evangelical",
                      "4.religion" = "other denomination",
                      "pope_francis" = "evaluation Pope Francis (H3)",
                      "0b.conf_chrch" = "no confidence at all (ref.)",
                      "1.conf_chrch"  = "little confidence",
                      "2.conf_chrch"  = "some confidence",
                      "3.conf_chrch"  = "a lot of confidence",
                      "1.female" = "female",
                      "leftright" = "left-right position",
                      "age_c" = "age",
                      "1.married" = "married",
                      "1b.race" = "black (ref.)",
                      "2.race" = "indigenous", 
                      "3.race" = "mestizo", 
                      "4.race" = "white",
                      "5.race" = "other",
                      "0b.relig_commit" = "not practicing at all (ref.)",
                      "1.relig_commit" = "not very practising",
                      "2.relig_commit" = "practising",
                      "3.relig_commit" = "very practising",
                      "1.education" = "illiterate",
                      "2.education" = "incomplete primary",
                      "3.education" = "complete primary",
                      "4.education" = "incomplete secondary",
                      "5b.education" = "complete secondary (ref.)",
                      "6.education" = "incomplete tertiary",
                      "7.education" = "complete tertiary")) %>%
  add_row(var = "Education", beta = 0) %>%
  add_row(var = "Trust in church (H2)", beta = 0) %>%
  add_row(var = "Religiosity", beta = 0) %>%
  add_row(var = "Ethnic group", beta = 0) %>%
  add_row(var = "Religious denomination (H1)", beta = 0) %>%
  add_row(var = " ", beta = 0) %>%
  add_row(var = "  ", beta = 0) %>%
  add_row(var = "   ", beta = 0) %>%
  mutate(var = factor(var, levels = c("complete tertiary",
                                      "incomplete tertiary",
                                      "complete secondary (ref.)",
                                      "incomplete secondary",
                                      "complete primary",
                                      "incomplete primary",
                                      "illiterate",
                                      "Education",
                                      "other",
                                      "white",
                                      "mestizo",
                                      "indigenous",
                                      "black (ref.)",
                                      "Ethnic group",
                                      "married",
                                      "female",
                                      "age",
                                      "left-right position",
                                      "very practising",
                                      "practising",
                                      "not very practising",
                                      "not practicing at all (ref.)",
                                      "Religiosity",
                                      "   ",
                                      "evaluation Pope Francis (H3)",
                                      "  ",
                                      "a lot of confidence",
                                      "some confidence",
                                      "little confidence",
                                      "no confidence at all (ref.)",
                                      "Trust in church (H2)",
                                      " ",
                                      "none",
                                      "other denomination",
                                      "evangelical",
                                      "roman catholic (ref.)",
                                      "Religious denomination (H1)")))

vec_vars1 <- c("roman\ncatholic",
              "evangelical")
vec_vars2 <- c("no\nconfidence\nat all",
              "little confidence",
              "some confidence",
              "a lot\nof\nconfidence")
vec_vars <- c(vec_vars1, vec_vars1, vec_vars2, vec_vars2, "very\nbad\n(0)", "very\ngood\n(10)", "very\nbad\n(0)", "very\ngood\n(10)") 
vec_sel <- c(rep(0,2), rep(1,2), rep(0,4), rep(1,4), rep(0,2), rep(1,2)) 
vec_add <- cbind(vec_vars, vec_sel)
  
df_prob <- read_csv("tables and figures/model_pred_prob.csv") %>%
  mutate(lower_95 = beta - 1.96 * se,
         upper_95 = beta + 1.96 * se, 
         lower_99 = beta - 2.58 * se,
         upper_99 = beta + 2.58 * se) %>%
  cbind(., vec_add) %>% 
  filter(!vec_vars %in% c("little confidence", "some confidence")) %>%
  mutate(vec_vars = factor(vec_vars, levels = c("roman\ncatholic",
                                                "evangelical", 
                                                "no\nconfidence\nat all", 
                                                "a lot\nof\nconfidence",
                                                "very\nbad\n(0)",
                                                "very\ngood\n(10)")))

df_prob_int <- read_csv("tables and figures/model_pred_prob_int.csv") %>%
  mutate(lower_95 = beta - 1.96 * se,
         upper_95 = beta + 1.96 * se) %>%
  mutate(pope_francis = seq(10,0)) 

df_catholics <- read_csv("data/catholics_LAC.csv")


# 2. generate figures ----
## 2.1 Figure 1 ----
fig1 <- ggplot(data = df_catholics, 
               aes(x = country, y = catholics)) +
  geom_bar(stat = "identity") +
  scale_y_continuous(name = "% Catholics\n",
                     breaks = c(0,20,40,60,80,100)) +
  scale_x_discrete(name = "",
                   guide = guide_axis(angle = 45)) +
  theme_classic(base_size = 40) +
  theme(plot.background = element_rect(fill = bg_colour),
        panel.grid.major.y = element_line(color = "grey30",
                                          size = 0.5,
                                          linetype = 2))

## 2.2 Figure 2 ----
fig2 <- ggplot(data = df, aes(beta, var)) + 
  geom_point(size = 4) +
  geom_point(data = filter(df, beta == 0), size = 8, colour = "white") +
  geom_linerange(aes(y = var, xmin = lower_95, xmax = upper_95)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(name = "",
                     limits = c(-0.83, 0.83), 
                     breaks = c(-0.75, 0, 0.75),
                     labels = c("\n-0.75\n(lower belief\nin climate change)",
                                "0",
                                "\n0.75\n(higher belief\nin climate change)")) +
  scale_y_discrete(name = "Coefficient estimates of model 3\n") +
  theme_classic(base_size = 40) +
  theme(plot.background = element_rect(fill = bg_colour))

## 2.3 Figure 3 ----
fig3a <- ggplot(data = filter(df_prob, vec_sel == 1), aes(vec_vars, beta)) + 
  geom_point(size = 6) +
  geom_linerange(aes(x = vec_vars, ymin = lower_95, ymax = upper_95),
                 size = 1.2) +
  geom_hline(yintercept = 0, linetype = "solid") +
  geom_hline(yintercept = 0.3753, linetype = "dashed") +
  geom_segment(x = 1, xend = 2, y = 0.25, yend = 0.25) +
  geom_segment(x = 3, xend = 4, y = 0.25, yend = 0.25) +
  geom_segment(x = 5, xend = 6, y = 0.25, yend = 0.25) +
  geom_segment(x = 1, xend = 1, y = 0.25, yend = 0.26) +
  geom_segment(x = 2, xend = 2, y = 0.25, yend = 0.26) +
  geom_segment(x = 3, xend = 3, y = 0.25, yend = 0.26) +
  geom_segment(x = 4, xend = 4, y = 0.25, yend = 0.26) +
  geom_segment(x = 5, xend = 5, y = 0.25, yend = 0.26) +
  geom_segment(x = 6, xend = 6, y = 0.25, yend = 0.26) +
  geom_text(x = 1.5, y = 0.22, label = "+3.2 %-points", size = 10) +
  geom_text(x = 3.5, y = 0.22, label = "-9.0 %-points", size = 10) +
  geom_text(x = 5.5, y = 0.22, label = "+6.2 %-points", size = 10) +
  geom_text(x = 1.5, y = 0.48, label = "Religious\ndenomination", size = 10) +
  geom_text(x = 3.5, y = 0.48, label = "Trust in\nchurch", size = 10) +
  geom_text(x = 5.5, y = 0.48, label = "Evaluation\nPope Francis", size = 10) +
  scale_y_continuous(name = "Predicted probability\nto strongly believe\n",
                     limits = c(0, 0.5)) +
  scale_x_discrete(name = "") +
  theme_classic(base_size = 60) +
  theme(plot.background = element_rect(fill = bg_colour))

fig3b <- ggplot(data = filter(df_prob, vec_sel == 0), aes(vec_vars, beta)) + 
  geom_point(size = 6) +
  geom_linerange(aes(x = vec_vars, ymin = lower_95, ymax = upper_95),
                 size = 1.2) +
  geom_hline(yintercept = 0, linetype = "solid") +
  geom_hline(yintercept = 0.1079, linetype = "dashed") +
  geom_segment(x = 1, xend = 2, y = 0.04, yend = 0.04) +
  geom_segment(x = 3, xend = 4, y = 0.04, yend = 0.04) +
  geom_segment(x = 5, xend = 6, y = 0.04, yend = 0.04) +
  geom_segment(x = 1, xend = 1, y = 0.04, yend = 0.05) +
  geom_segment(x = 2, xend = 2, y = 0.04, yend = 0.05) +
  geom_segment(x = 3, xend = 3, y = 0.04, yend = 0.05) +
  geom_segment(x = 4, xend = 4, y = 0.04, yend = 0.05) +
  geom_segment(x = 5, xend = 5, y = 0.04, yend = 0.05) +
  geom_segment(x = 6, xend = 6, y = 0.04, yend = 0.05) +
  geom_text(x = 1.5, y = 0.02, label = "-1.2 %-points", size = 10) +
  geom_text(x = 3.5, y = 0.02, label = "+3.1 %-points", size = 10) +
  geom_text(x = 5.5, y = 0.02, label = "-2.3 %-points", size = 10) +
  geom_text(x = 1.5, y = 0.48, label = "Religious\ndenomination", size = 10) +
  geom_text(x = 3.5, y = 0.48, label = "Trust in\nchurch", size = 10) +
  geom_text(x = 5.5, y = 0.48, label = "Evaluation\nPope Francis", size = 10) +
  scale_y_continuous(name = "Predicted probability\nto disbelieve\n",
                     limits = c(0, 0.5)) +
  scale_x_discrete(name = "") +
  theme_classic(base_size = 60) +
  theme(plot.background = element_rect(fill = bg_colour))

## 2.4 Figure 4 ----
fig4 <- ggplot(data = df_prob_int, aes(pope_francis, beta)) +
  geom_point(size = 6) +
  geom_linerange(aes(x = pope_francis, ymin = lower_95, ymax = upper_95), 
                 size = 1.5) +
  geom_hline(yintercept = 0, linetype = "solid") +
  geom_hline(yintercept = 0.3416, linetype = "dashed") +
  scale_y_continuous(name = "Predicted probability\nto strongly believe\n| a lot of confidence\n",
                     limits = c(0, 0.5)) +
  scale_x_continuous(name = "Evaluation\nPope Francis",
                     limits = c(0, 10.2),
                     breaks = c(seq(0, 10)),
                     labels = c("very\nbad (0)",
                                "1", "2", "3", "4", "5", "6", "7", "8", "9",
                                "very\ngood (10)")) +
  theme_classic(base_size = 60) +
  theme(plot.background = element_rect(fill = bg_colour))


# export figures
jpeg("tables and figures/fig1.jpeg", quality = 100, 
     height = 1200, 
     width = 1960,
     units = "px")
fig1
dev.off()  

jpeg("tables and figures/fig2.jpeg", quality = 100, 
     height = 1000, 
     width = 1600,
     units = "px")
fig2
dev.off()  

jpeg("tables and figures/fig3.jpeg", quality = 100, 
     height = 2000, 
     width = 1600,
     units = "px")
fig3 <- grid.arrange(fig3a,
                     fig3b,
                     ncol = 1, 
                     nrow = 2)
dev.off()  

jpeg("tables and figures/fig4.jpeg", quality = 100, 
     height = 1200, 
     width = 1960,
     units = "px")
fig4
dev.off()  
